Algorytm p-adaptacji
Celem algorytmu p adaptacji jest zwiększenie dokładności aproksymacji na siatkach obliczeniowych poprzez zwiększenie stopnia wielomianów rozpiętych na siatce obliczeniowej. Algorytm p adaptacji możliwy jest do zaimplementowania dopiero wtedy kiedy skonstruujemy takie funkcje bazowe rozpięte na elementach skończonych, dla których możliwe jest podnoszenie stopnia wielomianów. W przypadku tradycyjnych funkcji B-spline, opisywanych szerzej w tym podręczniku, możliwe jest jedynie globalne podniesienie stopnia funkcji B-spline na lokalnym patchu elementów skończonych.
Innymi słowy, załóżmy że mamy wektor węzłów opisujący bazę B-spline'ów drugiego stopnia o ciągłości C1 w kierunku osi \( x \) [0 0 0 1 2 3 4 4 4], oraz identyczny wektor węzłów opisujący bazę B-spline'ów drugiego stopnia w kierunku osi \( y \) [0 0 0 1 2 3 4 4 4].
Wektory te rozpinają dwie jednowymiarowe bazy funkcji B-spline w kierunku poszczególnych osi
\( B^x_{1;2}(x),...,B^x_{6;2}(x); B^y_{1;2}(y),...,B^y_{6;2}(y) \)
które poprzez iloczyn tensorowy definiują bazę dwuwymiarowych funkcji bazowych
\( \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=1,...,6;j=1,...,6 \)
Możliwe są następujące p adaptacje:
- podniesienie stopnia funkcji B-spline w kierunku osi \( x \) o 1, poprzez zmodyfikowanie wektora węzłów w kierunku osi \( x \) [0 0 0 0 1 2 3 4 4 4 4]. Dostaniemy wówczas bazy \( B^x_{1;3}(x),...,B^x_{7;3}(x); B^y_{1;2}(y),...,B^y_{6;2}(y) \) oraz dwuwymierowe funkcje \( \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=1,...,7;j=1,...,6 \)
- podniesienie stopnia funkcji B-spline w kierunku osi \( y \) o 1, poprzez zmodyfikowanie wektora węzłów w kierunku osi \( y \) [0 0 0 0 1 2 3 4 4 4 4]. Dostaniemy wówczas bazy \( B^x_{1;3}(x),...,B^x_{7;3}(x); B^y_{1;2}(y),...,B^y_{6;2}(y) \) oraz dwuwymierowe funkcje \( \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=1,...,6;j=1,...,7 \)
- podniesienie stopnia funkcji B-spline w kierunku obydwu osi, poprzez zmodyfikowanie wektora węzłów w obydwu kierunkach [0 0 0 0 1 2 3 4 4 4 4] oraz [0 0 0 0 1 2 3 4 4 4 4] o 1. Dostaniemy wówczas bazy \( B^x_{1;3}(x),...,B^x_{7;3}(x); B^y_{1;3}(y),...,B^y_{7;3}(y) \) oraz dwuwymierowe funkcje \( \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=1,...,7;j=1,...,7 \)
Możemy oczywiście podnosić stopnie funkcji bazowych w jednym lub w drugim kierunku o dowolny stopień (o ile nasz program komputerowy potrafi wygenerować B-spline wyższego stopnia bez problemu).
W przypadku funkcji B-spline nie ma możliwości mieszania stopni wielomianów bez umieszczania \( C^0 \) separatorów pomiędzy funkcjami B-spline. Możliwe jest natomiast mieszanie patchów (grup elementów) przedzielonych separatorami, na których wielomiany B-spline mają różną ciągłość.
Na przykład załóżmy że mamy wektor węzłów opisujący bazę B-spline'ów drugiego stopnia o ciągłości C1 na 2 elementach w kierunku osi \( x \) oraz wektor węzłów opisujący bazę B-spline'ów trzeciego stopnia o ciągłości C2 na kolejnych 2 elementach w kierunku osi \( x \) czyli [0 0 0 1 2 2 2] [2 2 2 2 3 4 4 4 4].
Wektory te rozpinają razem bazę
\( B^x_{1;2}(x),...,B^x_{6;2}(x); B^y_{1;2}(y),...,B^y_{6;2}(y) \)
Wektory te możemy "złożyć" z bazą B-spline'ów drugiego stopnia w kierunku osi \( y \) [0 0 0 1 2 3 4 4 4], czyli \( B^x_{1;2}(x),...,B^x_{6;2}(x); B^y_{1;2}(y),...,B^y_{6;2}(y) \).
Wektory te rozpinają bazę funkcji B-spline w kierunku osi \( x \) w następujący sposób
\( B^x_{1;2}(x),B^x_{2;2}(x),B^x_{3;2}(x),B^x_{4;2}(x),B^x_{5;3}(x),B^x_{6;3}(x),B^x_{7;3}(x),B^x_{8;3}(x),B^x_{9;3}(x) \).
Poprzez iloczyn tensorowy otrzymujemy bazę dwuwymiarowych funkcji bazowych
\( \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=1,...,4;j=1,...,6 \} \cup \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=5,...,9;j=1,...,6 \} \)
Zauważmy że grupy elementów sklejone w ten sposób nie dają ciągłej aproksymacji w punkcie 2. Z lewej strony mamy jeden wielomian osiągający maksymalną wartość 1 w punkcie 2, z prawej strony mamy drugi inny wielomian osiągający maksymalną wartość 1 w punkcie 2.
Para wektorów węzłów [0 0 0 1 2 2 2] [2 2 2 2 3 4 4 4 4] jest równoważna wektorowi [0 0 0 1 2 2 2 2 3 4 4 4 4]. Oczywiście algorytm interpretujący taki wektor węzłów musi zauważyć że począwszy od węzła 2 który powtórzony został 4 razy, zwiększa nam się stopień wielomianów o 1. Dodatkowo, ostatni wielomian w pierszej podgrupie elementów (wielomian \( B^x_{4;2}(x) \) należało by scalić z pierwszym wielomianem w drugiej podgrupie elementów \( B^x_{5;3}(x) \). Scalenie takie jest możliwe ponieważ pierwsze i ostatnie człony B-spline'ów są identycznymi wielomianami niezależnie od ich stopnia. Uzyskamy wtedy wielomian drugiego stopnia \( B^x_{4;2}(x) \) który rozpościera się na ostatni element w pierwszej podgrupie i na pierwszy element w drugiej podgrupie. Nasze bazy wyglądają więc teraz tak
\( B^x_{1;2}(x),B^x_{2;2}(x),B^x_{3;2}(x),B^x_{4;2}(x),B^x_{5;3}(x),B^x_{6;3}(x),B^x_{7;3}(x),B^x_{8;3}(x) \).
\( B^x_{1;2}(x),...,B^x_{6;2}(x); B^y_{1;2}(y),...,B^y_{6;2}(y) \)
\( \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=1,...,4;j=1,...,6 \} \cup \{ B^x_{i;2}(x)B^y_{j;2}(y)\}i=5,...,8;j=1,...,6 \} \)
- jednorodna p adaptacja wykonywana apropri (siatka początkowa, dokładność, maksymalna liczba kroków)
- aktualna siatka = siatka początkowa
- rozwiąż swój problem obliczeniowy na aktualnej siatce (na przykład projekcję bitmapy lub transport ciepła), obliczając aktualne rozwiązanie \( u \)
- pętla dla i = 1 do maksymalnej liczby kroków
- \( K_{maxerr }=0; K^x_{maxerr}=0; K^y_{maxerr}=0 \)
- dla każdego elementu \( K \) aktualnej siatki wykonaj następujące czynności
- oblicz bład aproksymacji na elemencie \( K \) obliczając seminormę gradientu rozwiązania na elemencie \( K_{err}=\|\nabla(B - u) \|=\int |\frac{\partial (B-u)}{\partial x}|^2+ |\frac{\partial (B-u)}{\partial y}|^2dxdy \), oraz normy z pochodnych kierunkowych \( K^x_{err}=\|\nabla_x(B - u) \|=\int |\frac{\partial (B-u)}{\partial x}|^2dxdy \) oraz \( K^y_{err}=\|\nabla_y(B - u) \|=\int \|\frac{\partial (B-u)}{\partial y}|^2dxdy \)
- jeśli \( K_{err}>K_{maxerr} \) wtedy \( K_{maxerr}=K_{err} \)
- jeśli \( K^x_{err}>K^x_{maxerr} \) wtedy \( K^x_{maxerr}=K_{err} \)
- jeśli \( K^y_{err}>K^y_{maxerr} \) wtedy \( K^y_{maxerr}=K_{err} \)
- koniec pętli po elementach
- jeśli \( K_{maxerr}< \) wymagana dokładność, to koniec
- pętla po elementach \( K \) siatki
- jeśli \( K_{err}>0.33*K_{maxerr} \) wówczas podnieś stopień aproksymacji wielomianowej w dwóch kierunkach, kontynuuj pętlę po elementach
- jeśli \( K^x_{err}>0.33*K^x_{maxerr} \) wówczas podnieś stopień aproksymacji wielomianowej w kierunku osi \( x \), kontynuuj pętlę po elementach
- jeśli \( K^x_{err}>0.33*K^y_{maxerr} \) wówczas podnieś stopień aproksymacji wielomianowej w kierunku osi \( y \), kontynuuj pętlę po elementach
- koniec pętli po elementach
- konieć pętli adaptacji
Jeśli dany element wymaga modyfikcji stopnia wielomianów w celu polepszenia jakości aproksymacji, możliwy jest jeden z trzech opisanych scenariuszy, podniesenie stopnia funkcji bazowych w kierunku osi \( x \), podniesienie stopnia funkcji bazowych w kierunku osi \( y \), lub podniesnie stopnia w obydwu kierunkach.
W algorytmie adaptacyjnym zazwyczaj nie wykonuje się adaptacji na wszystkich elementach siatki, ale jedynie na tych elementach na których błąd jest większy niż 33 procent błedu maksymalnego (liczonego po wszystkich elementach). Jest to oczywiście arbitralnie przyjęta wartość. W zaproponowanej wersji algorytmu próbujemy najpierw wykonać adaptację w obydwu kierunkach (podnieść stopień wielomianu o 1 w obydwu kierunkach), a następnie (jeśli nie jest to wskazane) próbujemy podnieść stopień wielomianów w jednym z dwóch kierunków. Jest to oczywiście pzykładowa wersja algorytmu, możliwe są różne modyfikacje. Najbardziej zaawansowaną wersją algorytmu adaptacyjnego jest algorytm hp adaptacji opisany w innym module.
Scalając ze sobą wiele patchów (grup elementów) na których rozpięte są różne wielomiany B-spline, uzyskamy siatkę obliczeniową na której możliwe będą lokalne modyfikacje stopnia funkcji B-spline, na przykład stosując opisany algorytm p-adaptacji.
Gdy patch elementowy składa się z jednego elementu, funkcje bazowe rozpięte na tym elemencie odpowiadają klasycznemu algorytmowi p-adaptacji pracującemu na wielomianach Lagrange'a. Dzieje się tak dlatego, iż wektory węzłów w których każdy element oddzielony jest \( C^0 \) separatem, tak jak ma to miejsce w klasycznych wielomianach Lagrange'a. Obie bazy generują wówczas identyczną przestrzeń wielomianów.
Matematyczne właściwości algorytmu p-adaptacji opisał po raz pierwszy prof. Ivo Babuśka z Uniwersytetu Teksańskiego w Austin
[1].